Towards a working density-functional theory for polymers: 
First-principles determination of the polyethylene crystal structure 
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Equilibrium polyethylene crystal structure, cohesive energy, and elastic constants are calculated by 
density-functional theory applied with a recently proposed density functional (vdW-DF) for general 
geometries [Phys. Rev. Lett. 92, 246401 (2004)] and with a pseudopotential-planewave scheme. The 
vdW-DF with its account for the long-ranged van der Waals interactions gives not only a stabilized 
crystal structure but also values of the calculated lattice parameters and elastic constants in quite 
good agreement with experimental data, giving promise for successful application to a wider range 
of polymers. 
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Understanding of crystalline solids is greatly enhanced 
by the periodicity of the atomic structure, which allows 
very detailed comparisons between theory and experi- 
ment. Macromolecular materials, such as polymers and 
other complex fluids, are less well understood, due to 
the challenging complex natures of their structures. Fea- 
tures of the polymer on the atomic and mesoscopic length 
scales are inter-dependent. For instance, atomic bonds 
are on the Angstrom scale, while diffusion processes in- 
volve whole chains extending some 100 A. In order to 
cover the full set of length scales by theory we need 
atomic-scale input for mesoscopic-scale force-scheme sim- 
ulations. Input forces may be determined routinely from 
first-principles density-functional theory (DFT) calcula- 
tions, the more accurate a functional the better, for a 
reasonably-sized unit cell. Among polymer crystals, with 
their typically very complex structure, the polyethylene 
(PE) crystal has a relatively simple structure [lj . A first- 
principles implementation of DFT 2], vdW-DF, with 
a general account for van der Waals (vdW) forces, is 
here shown to give results for crystal structure, cohe- 
sive energy, and elastic constants of the PE crystal in a 
very promising agreement with previous low-temperature 
measurements. 

Earlier attempts to study the equilibrium properties 
of polymer systems have relied on semiempirical inter- 
atomic potentials, including force-field methods [H, 0| 
and modified DFT calculations The fitting of in- 
ternuclear potentials in order to reproduce experimental 
equilibrium properties gives serious restrictions on the 
predictive power of these methods for nonequilibrium sit- 
uations and for molecules, for which experimental data 
are scarce or lacking. Further, these semiempirical meth- 
ods lack information to make a systematic improvement 
of the interatomic potentials possible. A physically mo- 
tivated first-principle description of sparse polymer crys- 
tals is highly desirable. When applied at perfect con- 
ditions, it might even render the possibility of develop- 
ing consistent and transferable interatomic potentials for 



both force schemes and hybrid methods. 

The ongoing DFT success story for ground-state elec- 
tron structure calculations of dense materials systems is 
driven by the relatively low computational cost together 
with the ability to describe very diverse systems. How- 
ever, the widely used local and semilocal implementa- 
tions do not include long-range nonlocal correlations that 
are essential for a proper description of the intermolecu- 
lar vdW interactions. These interactions are crucial for 
the stability of systems with regions of low electron den- 
sity, often encountered in biological and nanotechnologi- 
cal applications. The recently developed vdW-DF func- 
tional [2j accounts for the vdW interaction in a seamless 
way both at the equilibrium binding separation and at 
asymptotically large separations. It applies to general 
geometries and accounts for the nonlocal correlation via 
a functional that takes only the electron density as in- 
put. It has proved very promising in describing equilib- 
rium separations and binding energies for a range of sys- 
tems, including dimers of benzene rings 0, [f|, graphene 
sheets Q, polycyclic aromatic hydrocarbons (PAH's) [H, 
monosubstituted benzene molecules |9|, and combina- 
tions thereof, such as PAH's and phenol adsorbed on 
graphite @,[l(|. 

Simple polymer crystals constructed from long paral- 
lel polymer chains such as linear PE or polypropy- 
lene, represent an important class of sparse materials 
that are stabilized by the vdW interaction. At low tem- 
peratures, the PE crystal stabilizes in a simple base- 
centered orthorhombic crystal structure (Pna2i), Fig. [1] 
as established by X-ray and neutron scattering experi- 
Standard theory, first-principles 
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calculations with the semilocal generalized gradient ap- 
proximation (GGA) of DFT, predicts this structure to be 
unstable [B|, [l5| , calling for incorporation of nonlocal in- 
teractions. Recent vdW-DF calculations on parallel, well 
separated PE molecules [3, 17 1 show the vdW interac- 
tion in these systems to be non-negligible. To account 
for the nonlocal correlations at all separations, in partic- 
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FIG. 1: Schematics of the PE crystal structure in its base- 
centered orthorhombic unit cell with the polymers aligned in 
parallel. Lattice parameters a and b are indicated, while the 
parameter c, which is the polymer repetition length, is per- 
pendicular to the plane shown. Solid (broken) circles repre- 
sent CHb-units in the plane (a distance c/2 out of the plane). 
The angular orientation of the polymers is determined by the 
setting angle 9, defined as the angle between the a(> b) axis 
and the intra-polymer carbon plane. 



ular those relevant in the PE crystal, we here apply the 
general-geometry vdW-DF [2] . 

The structure of the isolated polymer has earlier been 
calculated in the GGA 17], and this structure is used 
as input for the crystal calculations (l8| . The internal 
geometrical structure of the isolated polymer chain thus 
fixes the repetition length c = 2.57 A of the PE crystal 
unit cell. 

Our vdW-DF calculations give binding-energy values 
at varying crystal-parameter settings, which are pre- 
sented as contour maps, and from which equilibrium val- 
ues for lattice parameters, cohesive energies, and elastic 
coefficients are extracted and compared with experimen- 
tal low-temperature data. Polymers might have several 
conformations with similar energies. Likewise, in a poly- 
mer crystal structure several local energy minima might 
exist. In particular, whereas the structure shown in Fig.[T] 
with a > b is experimentally observed, the structure that 
approximately corresponds to pulling one of the polymers 
half a unit length (c/2) out of the plane of the paper leads 
to essentially the same polymer chain packing and in turn 
to a similar configurational energy. In fact, the structure 
with a polymer pulled outward can also be described by 
the schematics in Fig. [1] when the requirement a > b 
is relaxed. Our first-principles method enables us to di- 
rectly assess the cohesive-energy contours in both (a, b, 
6) regions. 

The vdW-DF functional is derived and described in 
Ref.fl It divides the correlation-energy functional into 
local and nonlocal parts, 



E c « E, 



LDA 



E*r 



(1) 



where the first term is approximated by the local-density 
approximation (LDA) and the second term vanishes for a 
uniform system. The dispersive interactions give the sec- 




FIG. 2: The nonlocal correlation energy for pairs of PE. 
For the practical evaluation of E^ in the crystal the range of 
integration in Eq. (f2Jl is important. Simple pair interaction 
calculations give a fast estimate of the range needed, although 
the full calculation includes all density within the integration 
range. Shown here are the relative sizes (meV per polymer 
segment) of the polymers' contributions to E™ (as pairs with 
the polymer at the circle centers), evaluated at the separations 
(A) and relative angles shown. The contours show length- 
averaged electron densities. Polymers that are closer than 
9 A (thick dashed circle) to the center polymer contribute 
with an energy of 321 meV. This is 99.4% of the total nonlocal 
correlation energy E^ of all the polymers shown (323 meV). 



ond term a substantial nonlocality, which allows a sim- 
pler account of the polarization properties from which it 
originates. It is determined from the inverse dielectric 



function 



which is assumed to be dominated by a 



single pole, with a strength fixed by the f-sum rule. The 
functional form obtained approximately is 



dr'n(r)0(r, r')n(r') . 



(2) 



The nonlocal kernel 4> can be tabulated Q in terms of 
a dimensionless distance D = (qo + q Q )\r — r'|/2 and an 
asymmetry parameter 5 — (qo — q )/(qo + Qo), where qo 
is a local parameter that depends on the electron density 
and its gradient at position r. The quantity qo is related 
to the pole position in e _1 , which is determined by the 
requirement that this e should give the same appropriate 
semilocal exchange-correlation energy density component 
as that of an electron gas in a much better approximation. 

The nonlocal-energy integral has the electronic den- 
sity n(r) as input. We use a self-consistently determined 
GGA density with the revPBE exchange flavor [2I Hil. l2fj] . 
calculated for the full crystal unit cell using a plane- wave 
code 2l| with ultrasoft pseudopotentials [22j . The plane- 
wave basis set is truncated at 400 eV, and a (4 x 4 x 10) 
/c-point Monkhorst-Pack sampling of the Brillouin zone 
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is used for the the periodically repeated unit cell. A spa- 
tial sampling separation of 0.13 A between fast-Fourier- 
transform grid points is used. 

In Ref. the vdW-DF is described for and applied 
to finite molecules. The generalization to extended sys- 
tems is straightforward. For an adsorbate system (a fi- 
nite molecule plus an extended surface) this is described 
in Ref. 0- For a bulk system, like the PE crystal, E^ l [n] 
of the crystal unit cell must include the interactions from 
the surrounding cells. This is taken care of by letting the 
spatial integrals in Eq. @ (as energy per unit cell) ex- 
tend over the unit cell for r and "everywhere" for r'. In 
practice, the r' integral is carried out over a region in 
space that is sufficiently large, i.e. by adding more space 
we do not change E™ 1 at the desired level of accuracy. 
We find it sufficient to include a total of 5 unit cells in 
each of the a and b directions (Fig. []} and 9 unit cells in 
the c direction in the r' integral, by which is found to 
deviate less than one part per thousand from the results 
obtained by integrating over a substantially larger region 
of space. Figure[2]illustrates this by showing the pairwise 
contributions to the polymer interaction E~. The figure 
shows that the most important contributions come from 
the nearest and next-nearest neighbor polymers. 

The total energy of the crystal is thus 



-EvdW-DF = Eqqa — EqGA,c + ^LDA.c + E£ 



(3) 



with the revPBE flavor of GGA. The cohesive energy is 
obtained as the difference between the calculated total- 
energy values for the actual structure and for a refer- 
ence structure with widely separated PE polymers. The 
electron-density-grid spacing as well as the polymer po- 
sition relative to the grid are kept fixed for the vdW- 
DF reference calculation to ensure that any polymer 
self-interaction of the crystal calculation and most grid- 
related numerical errors are canceled by the reference cal- 
culation. The scheme applies also to other sparse-matter 
systems. 

The cohesive-energy contours of the PE-crystal, that is 
the loci of equal binding-energy values for configurations 
with varying values of the lattice parameters a, b and the 
setting angle 9 (Fig. [3]), are calculated with vdW-DF, and 
two nearly equally deep minima are found at two different 
setting-angle values (Table |T]). The configuration approx- 
imately sketched in Fig. Q] is the most stable of the two. 
However, the two configurations are very close in energy, 
differing by only 1 meV/CH 2 (1%), and an accuracy suf- 
ficiently high to resolve such a stability issue cannot be 
claimed for the vdW-DF method. 

The calculated values of equilibrium lattice parame- 
ters for the Fig. Q] structure are in close agreement with 
the experimental ones (Table HJ), the vdW-DF values 
(a , b ) = (7.30 A, 5.22 A) differing by at most 8% from 
the experimentally observed values. The calculated val- 
ues of the cohesive energy (not including zero-point fluc- 




FIG. 3: The cohesive-energy contour map calculated with 
vdW-DF, expressed as energy per CH2 group as a function 
of the lattice parameters a and b. The plot (but not the cal- 
culation) is restricted to show the energetics at setting angle 
e = 44° , which is the optimal angle found with vdW-DF (Ta- 
ble Hi. The unit-cell height is fixed to c = 2.57 A. The full 
(dashed) contour curves are separated by 1 (10) meV/CHb, 
with values given in meV/CH2 at selected contours. The equi- 
librium parameters (summarized in Table [Jl obtained from 
experiments and from a LDA-DFT study are shown for com- 
parison. 



tuations) and the elastic coefficients (Table are also in 
good accord with experiments. 

In comparison with other attempts to calculate the 
PE-crystal equilibrium data (Table fl}, the vdW-DF re- 
sults show a very promising agreement with experimental 
data. This is, indeed, in great contrast to the fact that 
the PE-crystal is unstable in GGA, according to earlier 
studies [f], [ID, [23| and confirmed here, and also to the 



23]. 



Calculations with 
27] 



clear overbinding of the LDA 
ad-hoc corrections for the dispersive interactions 
with the BLYP functional 5] and with the force-field 
method [3j give results in a fair agreement with exper- 
imental data. However, the former uses an empirically 
damped —R~ 6 internuclear potential [5j, and the latter 
has parameters fitted to experimental data at 4 K, such 
as lattice parameters [3J, providing agreement with ex- 
periment at this temperature by construction. As a con- 
trast, the vdW-DF calculations do not take any empirical 
input. 

The vdW-DF functional used here combines the 
correlation-energy functional E^ 1 (Eq. with an ex- 
change functional taken from the revPBE flavor of GGA. 
Although the latter closely represents the Hartree-Fock 
exchange at separations relevant for this work, it gives 
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TABLE I: Values of equilibrium lattice parameters ao and 60, setting angle #0, elastic constants On, C12, and C22, and cohesive 
energy i? co h per CH2 group, calculated by the vdW-DF method and compared with other theoretical and experimental data. 



Method 


Cn (GPa) 


C12 (GPa) 


622 (GPa) E co \ 


(eV/CH; 


) ao (A) 


00 (A) 


So (degr.) 


vdW-DF 


11.1 


6.3 


8.9 


0.106 J 


7.30 


5.22 


44 


vdW-DF 


8.8 


6.7 


12.8 


0.105 J 


5.29 


7.28 


51 


LDA 


At O 

41. 2 


11.3 


47.1 


n 1 1 7 

U.ll 


6.55 


4.45 


44.4 


LDA+BLYP+6-LJ 6 


14.1 


7.2 


11.8 


0.1F 


7.40 


4.90 


45 


Force field (4 K) c 
Neutrons (77K) d 


14.0 


7.9 


13.5 


0.081 


7.20 


4.80 


41.9 


11.5 














Experiment 6 








0.080 








X rays (77 K) f 










7.42 


4.96 


47.7 


X rays (77 K) 9 










7.39 


4.93 


45 


Neutrons (4 K) h 










7.12 


4.85 


41 



Experiment (213 K)' 8.4 4.2 



a Ref. 23. 6 Ref. 5, c Ref..3. "Ref. 24. e Ref. 25. J Ref. V2. 9 Ref. P3. ft Ref. 14. l Ref. 26. 

3 Value not corrected for zero-point motion, estimated to account for approximately 0.01 eV/CHh. 



a slight misrepresentation of the true exchange repul- 
sion [2I, [H, Q . An analysis of the PE-molecule density 
profile implies that this effect should be stronger in the 
b than in the a direction, which is indeed the case for 
our calculated results. Thus an improved exchange de- 
scription should give an even better agreement with the 
experimental geometry. 

In summary, the vdW-DF method successfully predicts 
structural, cohesive, and elasticity data for the important 
test case of the PE crystal. It is gratifying that the inclu- 
sion of the fully nonlocal vdW-DF correlation functional 
[Eq. p])]. which has no empirical input or fitted parame- 
ters, leads to a good agreement with experimental data. 
As similar conclusions are drawn for a variety of differ- 
ent carbon-based systems, like graphene, dimers of ben- 
zene, polyaromatic hydrocarbons, monosubstituted ben- 
zene molecules, and combinations thereof 0,0, [USES]; 
the vdW-DF functional is certainly very promising for 
general kinds of geometries and molecules, like those in 
macromolecular materials. 
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